function ok = A4a_elasticities(i_k)

KK = 4092;

cd 'D:\data_replication'

load('estimation/5_supply_side/main_data_supply.mat')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Import Data and Estimates for Product k %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k = index_price_coeff(i_k,1);    


data_k = data_summary{k};
numProdsTotal = size(data_k,1);
X_MPEC_k = X_MPEC_summary{k};
X_MPEC_opt = X_MPEC_opt_summary{k};
share_k = data_k(:,5);
price_k = data_k(:,7) + 1;
gdp_capita_k = data_k(:,13);
distance_k = data_k(:,14);
population_k = data_k(:,16);
population_log_k = log(population_k);
Y_k = Y_summary{k};
FE_k = FE_summary{k};
market_id_k = data_k(:,end);
T = market_id_k(end,1);
nn = size(Y_k,2);
alpha_constant = X_MPEC_opt(1,1);
alpha_distance = X_MPEC_opt(2,1);
alpha_population = X_MPEC_opt(3,1);
alpha_price = X_MPEC_opt(4,1);
alpha_1 = X_MPEC_opt(12,1);
alpha_FE = X_MPEC_opt((3*(K_A+1) + 1):(3*(K_A+1) + (K - K_A)), 1);

xi_k = X_MPEC_opt((3*(K_A+1) + (K - K_A) + 2 + 1):(3*(K_A+1) + (K - K_A) + 2 + numProdsTotal), 1) + 10;
xi_k = xi_k - alpha_constant - log(price_k)*alpha_price - alpha_population*population_log_k - alpha_distance*distance_k;
q_k = alpha_constant + xi_k + FE_k*alpha_FE;


prodsMarket = zeros(T,1);
prodsMarket_temp = data_k(:,17);
for t = 1:T
prodsMarket_temp2 = prodsMarket_temp(market_id_k==t);
prodsMarket(t,1) = prodsMarket_temp2(1,1);
end

marketStarts = zeros(T,1);
marketEnds = zeros(T,1);
marketStarts(1) = 1;
marketEnds(1) = prodsMarket(1);
for t=2:T
    marketStarts(t) = marketEnds(t-1) + 1;
    marketEnds(t) = marketStarts(t) + prodsMarket(t) - 1;
end

marketForProducts = zeros(numProdsTotal,1);
for t=1:T
marketForProducts(marketStarts(t):marketEnds(t)) = t;
end
numProdsTotal = size(data_k, 1);


delta_k =  alpha_constant + xi_k + log(price_k)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;

for t = 1:T   
Y_market = Y_k(t,:);
expmu(marketStarts(t):marketEnds(t),:) = exp(log(price_k(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn));
end
expmeanval = exp(delta_k);
oo = ones(1,nn);                  
sharesum = sparse(zeros(T,numProdsTotal));  % used to create denominators in logit predicted shares (i.e. sums numerators)
for t = 1:T
    sharesum(t,marketStarts(t):marketEnds(t)) = 1;
end
numer = (expmeanval*oo ).*expmu;        
sum1 = sharesum*numer;
sum11 = 1./(0+sum1);                                                       % Will not result in share_k = EstShare_true
denom1 = sum11(marketForProducts,:);    
simShare_true = numer.*denom1;              
EstShare_true = mean(simShare_true,2);  

Y_k(marketForProducts,:);

alpha_i = alpha_price + alpha_1 * Y_k(marketForProducts,:);

own_elast = sum([(-alpha_i + 1) .* simShare_true + alpha_i .* (simShare_true .^ 2)], 2) ./ EstShare_true;
own_elast  = own_elast  / 100;
own_elast_median = median(own_elast);


for i = 1:numProdsTotal   
    i
price_k_counter = price_k;
t = marketForProducts(i,1);
price_k_counter(marketStarts(t):marketEnds(t),:) = price_k_counter(marketStarts(t):marketEnds(t),:) * 1.00001;
price_k_counter(i,1) = price_k(i,1);

delta_k_counter =  alpha_constant + xi_k + log(price_k_counter)*alpha_price + alpha_population*population_log_k + alpha_distance*distance_k;

for t = 1:T   
Y_market = Y_k(t,:);
expmu_counter(marketStarts(t):marketEnds(t),:) = exp(log(price_k_counter(marketStarts(t):marketEnds(t),:))*alpha_1*Y_market + repmat(FE_k(marketStarts(t):marketEnds(t),:)*alpha_FE,1,nn));
end
t = marketForProducts(i,1);
expmeanval_counter = exp(delta_k_counter);
numer_counter = (expmeanval_counter*oo ).*expmu_counter;        
sum1_counter = sharesum*numer_counter;
sum11_counter = 1./(0+sum1_counter);                                                       % Will not result in share_k = EstShare_true
denom1_counter = sum11_counter(marketForProducts,:);    
simShare_counter = numer_counter.*denom1_counter;              
EstShare_counter = mean(simShare_counter,2);  

x_own_true = EstShare_true(i,1) ./ price_k(i,1);
x_own_counter = EstShare_counter(i,1) ./ price_k_counter(i,1);

x_others_true = sum((EstShare_true(marketStarts(t):marketEnds(t),:)) ./ price_k(marketStarts(t):marketEnds(t),:)) - EstShare_true(i,1) ./ price_k(i,1);
x_others_counter = sum((EstShare_counter(marketStarts(t):marketEnds(t),:)) ./ price_k_counter(marketStarts(t):marketEnds(t),:)) - EstShare_counter(i,1) ./ price_k_counter(i,1);

price_others_true = [sum(price_k(marketStarts(t):marketEnds(t),:)) - price_k(i,1)] / prodsMarket(t,1);
price_others_counter = [sum(price_k_counter(marketStarts(t):marketEnds(t),:)) - price_k_counter(i,1)] / prodsMarket(t,1);

eos_num = (x_own_counter/x_others_counter - x_own_true/x_others_true) / (x_own_true/x_others_true);
eos_denom = (price_others_counter/price_k_counter(i,1) - price_others_true/price_k(i,1)) / (price_others_true/price_k(i,1));

eos(i,1) = eos_num / eos_denom; 

end

eos_median = median(eos);


filename_matlab = sprintf('data/firm_size_distribution/results/elasticities/elast_%i_k.mat', i_k);  
save(filename_matlab, 'i_k', 'eos_median', 'own_elast_median');  

ok = 1

end

